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The influence of Rashba spin-orbit interaction on the spin dynamics of a topologically disordered 
hopping system is studied in this paper. This is a significant generalization of a previous investiga- 
tion, where an ordered (polaronic) hopping system has been considered instead. It is found, that 
in the limit, where the Rashba length is large compared to the typical hopping length, the spin 
dynamics of a disordered system can still be described by the expressions derived for an ordered 
system, under the provision that one takes into account the frequency dependence of the diffusion 
constant and the mobility (which are determined by charge transport and are independent of spin). 
With these results we are able to make explicit the influence of disorder on spin related quantities 
as, e.g., the spin life-time in hopping systems. 
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I. INTRODUCTION 

Much attention has recently been devoted to the prob- 
lem of utilizing electron spin in semiconductor electron- 
ics. An overview of this evolving subject of spintronics 
is given in Ref. Q. One central aim of these efforts is 
the (preferably) electrical generation, manipulation, and 
detection of spin currents. 

While a large part of spintronics related literature 
is concerned with itinerant electron systems, the corre- 
sponding behavior of localized electrons is worth serious 
consideration, too. For instance, some proposals for solid 
state quantum computing are based on manipulating the 
electronic or nuclear spin state of impurities. 2 - 3 - 4 - 5 

One possibility of affecting the spin behavior through 
electrical means is by utilizing two-dimensional struc- 
tures (semi-conductor interfaces or heterostructures) 
showing Rashba spin-orbit interaction (SOI). 6 This paper 
considers a model, which consists of a two-dimensional 
system of localized electronic states. Charge trans- 
port is possible through thermal activation (hopping 
transport) i The spin dynamics is assumed to be deter- 
mined by Rashba SOI. 

In a former publication^ we derived microscopic and 
macroscopic transport equations for charge and spin of 
an ordered (polaronic) hopping system with Rashba SOI. 
Here, we generalize this theory to the case of spatial 
(i.e. not energetic) disorder. A suitable reference system 
would be nearest-neighbor hopping between donors in 
a semi-conductor or between Anderson localized states. 
The central approximation concerning spin behavior in 
hopping systems, which applies to this investigation, con- 
sists in the assumption, that the spin degree of freedom 
is solely influenced by Rashba SOI, i.e., there must be 
no spin dependent scattering (only s-wave scattering, no 
magnetic scatterers). Furthermore, the Rashba length, 
i.e. the length scale over which the spin rotation occurs, 
must be large compared to a typical hopping length R t . 

In Sec. [H] the (disorder averaged) transport equations 
for spin are derived for a disordered hopping system 
and their relation to the charge transport equations is 



studied. The (averaged or macroscopic) evolution equa- 
tions governing charge transport (drift-diffusion equa- 
tions) keep their functional shape in the transition from 
an ordered to a disordered system, the only effect being, 
that the diffusion constant and the mobility become fre- 
quency dependent (or, expressed differently: the macro- 
scopic evolution equations for the disordered system have 
memory, even if the microscopic equations did not). 

We show, that the equations for spin dynamics due 
to Rashba SOI likewise keep their functional shape in 
the presence of disorder, provided the hopping length is 
sufficiently short. Furthermore, the only differences to 
the ordered case are again a frequency dependence of the 
diffusion constant and the mobility of the charge. Thus, 
provided the length scale for spin inversion due to Rashba 
SOI (which will be denoted as the "Rashba length" in the 
following) is large in comparison with the typical hopping 
length, the influence of disorder on the spin dynamics is 
entirely determined by the corresponding change of the 
charge transport coefficients. 

Two different situations are considered in detail: The 
temporal evolution of an initial spin polarization is in- 
vestigated (Sec. lIII|l and the stationary state of a system 
with spin polarized boundary conditions is determined 
(Sec. lIVf) . In both cases, the analytical treatment is com- 
pared with numerical simulations. 

It is found, that the spin life-time strongly increases 
with increasing disorder. On the other hand, the spa- 
tial behavior of the diffusing spin in the stationary case 
is practically unaffected by the introduction of disorder. 
Additionally, the finding, previously made for an ordered 
Rashba hopping system, that there is a critical electrical 
field, above which the total spin polarization has an oscil- 
latory component, is also stable against the introduction 
of topological disorder. 



Scction lVIl gives a summary of the results and discusses 
the experimental relevance. The two appendices present 
details of the derivation of the macroscopic spin evolution 
equation and its long wave-length limit. 
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II. THE TRANSPORT EQUATIONS 

The equations which govern the evolution of charge 
and spin density of a hopping system with Rashba 
SOI have been derived in Refs. |8||3- These microscopic 
transport equations are rate equations for the occupa- 
tion numbers and the spin expectation values of the 
localized sites. ' It is found, that the transition matrix 
elements between different localized sites obtains the 
spin-dependent factor exp (— ier ■ (K x Rm'm)) due to 
Rashba SOI. Here, <x denotes the vector of Pauli spin ma- 
trices, R, m ' m = R m ' — R m is the distance vector between 
sites m! and m, and K is a vector normal to the two- 
dimensional plane, the length of which is proportional to 
the Rashba coupling strength. Its inverse 1/K = 1/|K| 
has the dimension of a length and will in the following 
be called Rashba length. 

The derivation of the rate equations is subject to the 
following restrictions: (i) single particle approximation 
(i.e. no electron-electron interaction), (ii) Markovian rate 
equations (i.e. the relevant times are beyond the quantum 
kinetic regime), and (iii) the Rashba length is large com- 
pared with typical hopping length scales, KR t « 1. In 
this paper, we furthermore restrict ourselves to: (iv) two 
site hopping probabilities. This gives the leading contri- 
bution to hopping transport, the small parameter being 
the ratio between the resonance integral and the ener- 
getic width of the set of electron states. Ohmic charge 
transport is adequately described by this approximation, 
but higher order effects, as e.g. the Hall effect, are ne- 
glected. 

The corresponding rate equations read* 



dt 

mi 

for the occupation numbers and 
d 



(1) 



dt 



Prn = {^mim ' P mi W mim - p m W mmi j (2) 



for the spin orientation. Here, p m is the occupation num- 
ber of site m, p m is the expectation value of the spin op- 
erator on this site, W mim is the transition rate between 
sites mi and m, and D mim is a 3 x 3-matrix, which de- 
scribes a rotation of the spin during the transition. An 
explicit expression of this matrix is given in Eq. (|B1|I . 
Note, that the two sets of equations are decoupled in the 
chosen approximation (due to the restriction to two site 
hopping probabilities). This means, that the charge (or 
particle) transport has no influence on the spin transport 
and vice versa. Thus, such an interesting phenomenon as 
the spin Hall effect ji2*ii*i2i which already has been stud- 
ied for ordered hopping systems^ and which has sparked 
controversial discussions recentlyj 13 i 14 i 15 i 16 i 17 lies outside 
the scope of the theory presented here. 

Reference utilized these rate equations in order to 
study the behavior of an ordered (polaronic) hopping sys- 
tem. Here, a topologically disordered hopping system is 



the subject of investigation. Thus, a disorder average has 
to be performed, in order to obtain macroscopic trans- 
port equations. 

AppendixlAleives the technical details of the derivation 
of the macroscopic transport equations subject to the 
restrictions given above. It proceeds along the lines of 
a "generic" averaging procedure, following the approach 
used in Ref. UJI 

Since the disorder averaged system is homogeneous, it 
is convenient to work in wave-vector space. The disor- 
der averaged rate equations in Fourier-Laplace space are 
obtained as 



[s + W{s\Q)-W{s\q)]p(s\q)=p Q {q) 



(3) 



and 



s + W(s\0)- / d z qi W{s\q - qrMqi) 



p( s lq) 
p (q)- (4) 



Here, s is the Laplace variable (related to time) and q is 
the (in-plane, i.e. two-dimensional) wave vector. 

In the long wave-length limit the transition rates of an 
isotropic system can be expanded in the form 



W(s\q) = W(s\0) - D(s)q 2 - i/i(s)q ■ E, 



(5) 



where the (generally frequency dependent) in-plane dif- 
fusion constant D and mobility p have been introduced. 
The quantity E is the in-plane electric field. In a system 
of reduced symmetry, D and p, are tensors, but this case is 
not considered here. Note, that the asymmetry between 
in-plane and normal-to-planc diffusion and drift does not 
show here, since the wave vectors q are restricted to two 
dimensions. 

As shown in App. [B]the support of D(q) lies within 
q < 2K. Thus, provided K is sufficiently small, the long 
wave-length limit is also applicable to the integrand in 
Eq. (0J. The condition KR t <C 1 already assures this, 
since wave-vectors corresponding to the (inverse of the) 
typical hopping length arc within the long wave-length 
limit. This can safely be assumed, because already for 
charge transport, we are not concerned with effects due 
to higher than the second spatial derivative of p(r) in the 
corresponding evolution equation. 

Using expression JSJ, the integral in Eq. (0} can be 
transformed into the 3 x 3-matrix 



(W(0) - Dq 2 - ipq ■ E)D(r) 

- (2i£>q - fiE) ■ V o D(r) + DAD(r) 



(6) 



where the s-dependence has been dropped from the no- 
tation. The symbol o denotes the dyadic product. Thus, 
only the first three moments of the rotation matrix D(q) 
(expressed here as spatial derivatives at the origin of 
space) enter the calculations in the long wave-length 
limit. 
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Finally, inserting the rotation matrix D as given by 
Eq. (|_B1|) one obtains the rate equations in the long wave- 
length limit (see Add. IBlfor details of the calculation, I3 
is the 3x3 identity matrix) 



and 



D{s)q 2 + ip(s)q • E p(s\q) = p (q) 



D(s)q 2 + ip,(s)q • E + AD{s)K 2 ) / 3 



(7) 



+ 4D(s)KoK ■ p(s\q) 
+ [K x (4iD(«)q - 2/i(s)E) ] x p( a |q) - p (q). (8) 

These evolution equations agree exactly with the cor- 



responding equations of an ordered hopping system, de- 
rived in Ref. 0| except that here the diffusion constant 
D and the mobility p are frequency dependent. This 
frequency dependence is entirely determined by charge 
transport and does not depend on SOI within the ap- 
proximations considered here. Thus, one finds that un- 
der the condition KR t <C 1, the spin dynamics is affected 
by disorder only through the change of D and p. Note, 
that the above derivation is independent of the approx- 
imations used to derive W(s|q) (the procedure used to 
obtain concrete expressions for D(s) and p(s)). In par- 
ticular, it is not restricted to continuous time random 
walk, the procedure further considered in Ref. even 
though we followed the general approach used there. 

For reference, the transport equations and 
transformed to real space- and time-coordinates read 



dt 



p(t\r) 



(Mi 



D(t - *i)Ap(ii|r) - p(t - <i)E • Vp(ii|r) 



(9) 



J dh{D(t - t x )Ap(ti|r) - p(t - ix)E • Vp(ii|r) 

- AK 2 D{t - h)p(h\r) - AK 2 D(t - tx)e z p z {tx\v) 
- 4£>(i - h)K [V Pz (h\r) - e z V ■ p(h\r)} - 2p(t - h)K [e z E ■ p(h\r) - Ep z {h\r)] }. (10) 

I 



III. TOTAL SPIN EVOLUTION 

In this section we consider the temporal evolution of 
the total spin magnetization. In this case we have to set 
q = in Eqs. and 0. Then, Eq. immediately 
yields particle number conservation. Writing the corre- 
sponding equation for spin polarization in matrix form, 
while fixing the coordinate system such that E || e x and 
K II e,, one obtains 



s + ADK 2 -2pKE ' 

s + 4DK 2 
2pKE s + 8DK 2 



p(s)=p . (11) 



One can see, that, in zero electric field E = 0, the spin 
components decay with two different time constants (tak- 
ing D for the moment as frequency independent): The 
^-component with n = 1/8DK 2 , and the in-plane com- 
ponents with twice this value T2 = l/ADK 2 , i.e. the spin 
component perpendicular to the plane decays two times 
faster than the in-plane components, a fact which has 
also been found for ordered hopping systems^ and itiner- 
ant electrons Note, in particular, that the spin life- 
time is inversely proportional to the diffusion constant. 
Since D substantially decreases with increasing disorder 
in hopping systems, the spin life-time will strongly in- 



crease with increasing disorder. 

The frequency dependence of D complicates the mat- 
ter, but its overall effect is to further increase the de- 
cay time constant for later times (see the discussion in 
Sec. 01. Thus, the spin decay slows down further in the 
progress of time. 

Even in a finite electric field, the in-plane component 
of p perpendicular to the field follows a decay law [p y in 
Eq. Qllf)]. whereas the other two components are coupled 
and have the solution 



with the matrix 



1 



det(n) 



n- 



PxO 
PzO 



n 



s + 8DK 2 
-2pEK 



2pEK 
s + 4DK 2 



(12) 



(13) 



When D and p do not depend on frequency, this corre- 
sponds to a sum of exponential functions in time, so that 
the spin components arc cither hyperbolic or trigonomet- 
ric functions of time, times an exponential decay factor 
(see Ref. @). Here, in the disordered case, D and p are 
frequency dependent, and the transformation of the solu- 
tion of Eq. I|llf) into i-space is impossible without choos- 
ing a specific model for the frequency dependence D(s) 
and p(s) beforehand. 
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For large times, D and p approach their respective DC- 
values, i.e., the asymptotic behavior of p for large times 
is easily obtained. On the other hand, due to the expo- 
nential decay, the large time behavior is only relevant, 
if the time constant for DC-behavior is smaller than the 
time constant for spin decay 1/ADK 2 . 

The solution for initial condition p — e z with time- 
independent D and p is given for future reference 



Px(t) 



— 6DKt 



:Sm(2DK 2 \/e 2 - It) (14) 



-6DK t 



cos(2DK 2 \/e 2 - It) 



1 



sin(2L>^Ve 2 - 1*) 



(15) 



The dimensionless electric field is e = pE/(DK). This 
quantity is time-independent if the Einstein relation be- 
tween D and p is valid, even when D and p themselves 
depend on time. In this case, e furthermore does not 
depend on the disorder, which also strongly affects the 
quantities D and p themselves. 

The solution here is written with trigonometric func- 
tions, which is the appropriate form for e > 1. In the 
case e < 1, hyperbolic functions have to be used instead. 
Thus, the dimensionless electric field e discriminates be- 
tween two different behaviors of the total spin polariza- 
tion: exponential decay in small electric fields and an 
additional oscillation in large electric fields. The occur- 
rence of one or the other regime can be tuned by varying 
(one or several of) the electric field, the temperature, and 
the Rashba length. 

Note, that the conclusions of this section remain the 
same, whether one considers a single polarized spin ini- 
tially placed at the origin, or a homogeneously polarized 
system. 



IV. STATIONARY STATE WITH BOUNDARY 
CONDITIONS 

As a second case, we determine the stationary state of 
a system with an in-plane electric field and spin injection. 
Taking, as before, E oc e x , and boundaries parallel to the 
y-axis, the charge and spin densities can only depend on 
the ^-coordinate. Denoting the derivative with respect 
to a; by a prime, one obtains 



= D oP "(x) - p Ep'(x) 



(16) 



= D Q p"(x)-poEp'(x)~AD Q K 2 p(x)-AD Q K 2 e z p z (x) 
-AD Q K[ exP ' z {x)-e z p' x {x)] 

- 2p a KE [e z p x (x) - e x p z (x)} . (17) 

Here, Do = / °° dtD(t) and /io = / °° dtp(t) are the DC- 
values of the corresponding quantities. 

Specifically, if one considers a half-plane, and takes the 
boundary conditions at x = and x — oo to be p(0) = 



p(oo) = po, p(0) = e z , and p(oo) 
reads p(x) — po, 



0, the solution then 



p x (x) = e x / x A x sm(x/A), 
p y (x) = 0, 

Pz (x) = e- x/x [cos(x/A)- A z sm{x/A)}, (18) 

where the dimensionless electric field e = poE / (DqK) 
assumes the DC- value and is the parameter determining 

the quantities lu± = 



512, the 



amplitudes A x 



-32 



and A, 



the decay length A = 2(uj + — e)^ 1 /K and the oscillation 
length A — 2/ {lo-K). Note, that the disorder enters only 
through the ratio po/D a . Thus, provided the Einstein 
relation between Dq and /zo is valid, the spatial behavior 
of p{x) does not depend on the disorder, the relevant 
length scale being determined mainly by K, since the 
(dimensionless) quantity lu- only very weakly depends 
on its sole parameter e (specifically, 3.95 < < 4). 



V. NUMERICAL SIMULATION 

Two kinds of numerical simulations are performed. 
First, the temporal evolution of a single spin is deter- 
mined (calculation A, corresponding to Sec. IIII|) by solv- 
ing the differential equations. This yields, after ensemble 
averaging, the temporal behavior of an initially local- 
ized z-spin. Secondly, the system is subjected to bound- 
ary conditions, which correspond to spin injection, and 
the spatial behavior of spin polarization in the stationary 
state is determined (calculation B, corresponding to Sec. 
II V|) . This is done by solving a linear system of equations. 
Again, a disorder average must then be performed. The 
disorder averages are performed over about 50 or 1000 
realizations for calculation A or B, respectively. In each 
case, the statistical error is of the order of the symbol 
size in the figures. The number of sites in the calcula- 
tions varies between about 1600 and 5000. It has been 
taken care, that the results are not affected by boundary 
effects. 

The transition rates for the simulation are taken to be 



W„ 



v e 



-aR„ 



(19) 



where vq is the attempt frequency (subsuming, e.g., the 
temperature dependence of the rates) and determines the 
numerical time scale. The parameter a is the inverse 
localization length of a single impurity state. Denoting 
by JV the (areal) density of localized states, the quantity 
aAf^ 1 / 2 provides a measure of the disorder. 

Calculations A and B complement each other. Calcu- 
lation A allows to study the temporal behavior, but the 
numerical complexity forbids to investigate systems with 
large disorder (large fx/V -1 / 2 ). Large disorder implies a 
large dispersion (typically many orders of magnitude) of 
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relevant hopping times. Since the numerical time step 
for a differential equation solver has to be related to the 
smallest relevant time scale, but the long time physics is 
related to the largest relevant time scale, the study of sys- 
tems with large disorder leads to an enormous increase in 
the number of time steps to be calculated, if one is inter- 
ested in the long time behavior. Furthermore, for larger 
disorder, the typical behavior can only be found in sys- 
tems, which are spatially more extended. This increases 
the numerical expense even further. Thus, only relatively 
small values of aA/" -1 / 2 and/or tu can be investigated by 
a numerical integration of the rate equations. 

On the other hand, calculation B immediately gives the 
t — ► oo behavior with a reasonable effort even for systems 
with large disorder. However, the temporal behavior is 
not accessible. 

To be specific, in calculation A we take the initial con- 
dition p = e z . That means, we consider a system with 
an initial spin polarization in the z-direction. The three 
components of the calculated total spin polarization are 
shown in Fig. E] for the disorder parameter uN^ 1 ! 2 = 3. 
Figures [3 and |3 give the z-component for several disorder 
parameters. 

For a comparison with the analytical expressions, the 
time-dependent diffusivity and mobility have to be de- 
termined. In method A, where the time evolution of 
each system is known, effective quantities are deter- 
mined by D cS (t) = <r 2 )/(4i) = (£ m p m (t)R 2 J / (U) 
and pes(t) = (x)/(Et) = (^ m Pm(*)R-m • e x )/(Et) from 
charge diffusion and drift. The brackets (• • • ) denote the 
disorder average. Method B poses a problem here, since 
the diffusion co-efficient cannot be directly measured in 
this case. Fortunately only the quotient p/D, which 
could be obtained using the Einstein relation, enters the 
analytic expressions. 

The relation between the effective quantity D c g(t) and 
those occurring in the analytical expressions [D(t), the 
inverse Laplace transform of D(s)] is unfortunately not 
simple, but rather involves differentiation operations: 
D(t) — -^2{tD cf [(t)]. A corresponding expression con- 
nects /Lt e ff and \x. Numerical differentiation of noisy data 
(the disorder averaged D cS ) is highly unstable, so that 
the direct comparison of the numerical results and the 
analytic expressions (but using the numerically obtained 
D and fi) is not reasonable. The situation were much im- 
proved, if a general model for the behavior of D(t) for the 
system under consideration were available, but no such 
model is available for the range of parameters considered 
here. 

Thus, we take a different approach. The asymptotic 
values of D e f[(t) and J Q dt\D(ti) are the same for very 
small (tfo < 1) or for very large times (DC-behavior). In 
between, for moderate times, the detailed behavior dif- 
fers, but, since they change monotonically over several 
orders of magnitude^ (see Fig. 0J, we expect the Marko- 
vian limit of the transport equations to be valid in the 
following sense. Because D(t) and p(t) decrease over 
many orders of magnitude with increasing t, the main 



contribution to the time integrals in Eqs. © and i|10|) 
comes from small times t < 1/vq. The charge and spin 
density is taken out of the integral (as in the derivation 
of a Markovian equation), the remaining time integrals 

replaced by J dtD(t) w D e g(t), and for p accordingly. 
The resulting equations are solved analytically, but using 
the numerically obtained (from charge transport) quanti- 
ties D c f[(t) and p c g(t). The analytic solutions calculated 
in this way are displayed in Figs. Q to for compari- 
son. One can see, that the agreement with the numerical 
simulation is reasonable. 

This discussion can be somewhat improved, without 
the need to introduce a specific model for D(t). Figure 
0] shows, that approximately 

{Ax>, {vot < 1); 

DoaM- K , (vot > 1); (20) 
D , M>1), 

with constants (infinite frequency) and Do (DC 

value). (The third case cannot really be deduced from 
Fig. 01 but it is justified for any system which does not 
have an infinite memory.) The exponent is about k as 1/4 
for aAi- 1/2 = 3 and k « 1/2 for aA/"~ 1/2 = 7. Small and 
large times [the first and third case of Eq. (|20[l] thus per- 
mit the approximation f dtiD(ti) « D c g(t). On the 
other hand, for medium times [the second case of Eq. 
(|20|l . which corresponds to the relevant time scale of the 
numerical calculations] , one obtains 

/ dhDih) W (1 - K)Deff(t). (21) 

Jo 

Thus, the effective D for the analytic expression has to be 
slightly reduced (by the factor 1 — k). This has the effect 
[see Eqs. (|14|) and (|15fl ] of decreasing the decay rate and 
also the frequency, in effect moving the zero crossings 
to larger times. Calculations show, that this model of 
D c g(t) is too crude to give a better agreement with the 
numerics compared with the simple model advertised in 
the previous paragraph. The general trend of diminishing 
the decay rate and the frequency, however, is in accord 
with the numerical evidence (see Fig.^l. 

Note, that the oscillatory behavior of the total spin po- 
larization for e > 1, which has been predicted for ordered 
hopping systems^ survives the introduction of disorder. 
Furthermore, Figs.|2]and|3]show, that increasing disorder 
(increasing <x/V -1 / 2 ) leads to an increase of the spin life- 
time. This can be readily explained by noting the strong 
decrease of the diffusion constant in the transport equa- 
tions, a conclusion, which is confirmed by the agreement 
between numerical simulation and analytical calculation. 

The results of calculation B (spin injection with bound- 
ary conditions) are shown in Fig. |SJ The only parameter 
of the analytical expressions is the ratio Do/ po which is 
assumed to take the value given by the Einstein relation 
(kT/e). The agreement to the analytical expressions i|18|) 
shown in Fig. [S] is convincing. Note, that in this case, 
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no parameter had to be obtained from charge transport. 
Instead, the temperature used in the numerical simula- 
tion is used in the analytical calculation, too. Figure El 
shows only the results for the single disorder parameter 
aN~ x l 2 = 7, but the good agreement between numerics 
and Eq. IjlSjl holds also for the other investigated param- 
eters (disorder 3 < aN^ 1 / 2 < 7 in small or medium fields 
< e< 1). 

A deviation is only found in very large electrical fields 
(e ^> 1, see Fig. [SJ. The spatial spin coherence in cal- 
culation B diminishes with increasing <x/V~ 1//2 in large 
electrical fields. The numerical data can still be fitted 
to the analytical model, but with a larger value of e as 
given by the Einstein relation (not shown). This indi- 
cates, that the ratio Do/fJ-o increases beyond the equilib- 
rium value kT/e in large fields, but otherwise, the spin 
behavior is still described by Eq. I|18fl . This is indicative 
of a transition from isotropic to directed percolation in a 
growing electric field. This transition leaves the diffusion 
constant nearly unaffected, while the mobility decreases. 
This scenario can explain the observed behavior qualita- 
tively, but further research is necessary to illuminate the 
true physical basis of the observed behavior and exclude 
other possible explanations. 



VI. DISCUSSION 

We have derived macroscopic spin transport equations 
for spatially disordered hopping systems with Rashba 
SOI. It is found, that the introduction of disorder leaves 
the vectorial structure of these equations intact, the only 
effect being that diffusion constant and mobility become 
frequency dependent (or, expressed differently, that the 
transport equations obtain memory) . This frequency de- 
pendence is already determined by the charge transport 
behavior and does not depend on the specifics of spin 
transport. The derivation of this relation between or- 
dered and disordered hopping spin transport is subject 
to the approximation, that the Rashba length (the length 
scale of spin precession) is large against a typical hopping 
length. Furthermore, only two-site hopping processes can 
be dealt with in the present treatment, thus excluding, 
e.g., the discussion of a possible spin Hall effect. 

Other effects, previously predicted for ordered spin 
hopping, also occur for disordered spin hopping. For in- 
stance, the spin decay is exponential in small electric 
fields, whereas it obtains an oscillatory component in 
large electric fields. Thus, there is a finite critical field, 
dividing both regimes, also in the disordered case. It is 
also conceivable that the behavior "oscillates" between 
exponential and vibrational behavior. But this can be 
excluded, so far as D{t) and /i(i) change monotonically 
with time. 

In comparison to the ordered case, the spin life-time 
is shown to be strongly increased by the introduction of 
disorder. This is explained by the significantly reduced 
diffusion constant, which enters the life-time reciprocally 



and which sharply decreases with increasing disorder. 

For the stationary state of spins injected through a 
boundary into a topologically disorder hopping system, 
it is found, that the spatial behavior of the spin polar- 
ization is largely unaffected by varying the disorder. An 
exception to this insensitivity is the systematic reduc- 
tion of spin polarization with increasing disorder in large 
in-plane electric fields e> 1. But even here, the numeri- 
cal data indicate, that the phenomenological description 
of spin transport, which is derived in this paper, is still 
valid. The deviation appears to be due to the invalidity 
of the Einstein relation between diffusion co-efficient and 
mobility in strong electric fields for the considered sys- 
tem, and is not a consequence of the peculiarities of spin 
transport. 

A comment concerning the required smallness of the 
Rashba interaction is in order here. For the value K = 
Vf7 /10 (the largest value used in the simulations) the 
spin reverses over a distance of the order of about ten 
times the mean distance between impurities. The typ- 
ical hopping length also is several mean distances long. 
Thus, the central condition KR t < 1 is hardly valid. 
This (large) value of K has been taken for numerical 
reasons, but even so, the agreement with theory is very 
good. This shows, that the phenomenological descrip- 
tion of spin transport derived above is quite robust even 
for (relatively) large values of the Rashba interaction 
strength. Smaller (more realistic) values of the Rashba 
coupling K render the approximation all the more valid. 

Finally, a few words concerning the experimental rele- 
vance of this work. Usually given values for the Rashba 
SOI strength (1CT 9 eV cm)2ii2£ lead to a Rashba length 
of the order of l/K ss 100 A, if one assumes that the effec- 
tive electron mass is equal to its elementary massif This 
length scale is in the range of typical hopping lengths. 
Thus, depending on the system under consideration, the 
theory presented here is a candidate for a description 
of the physical behavior. For the value of the Rashba 
coupling given above and at a temperature of 10 K, the 
critical electrical field e — 1 , which differentiates between 
exponential and oscillatory behavior of the total spin po- 
larization, corresponds to a field of about 10 3 V/cm in 
natural units. 

To the authors knowledge, the value of the Rashba 
coupling strength in hopping systems has not been in- 
vestigated up to now. It may well be, that this value 
is smaller than for itinerant electrons, in which case the 
estimate given above for the Rashba length (100 A) must 
be increased. Then, it is even easier to comply with the 
condition KR t C 1. A further advantage of this case 
would be, that a modification of the Rashba coupling 
strength by an external gate electrode^ would require 
smaller applied voltages. 
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APPENDIX A: THE MACROSCOPIC 
TRANSPORT EQUATIONS FOR SPIN DENSITY 

The aim is to derive macroscopic transport equations 
starting from the microscopic equations l[T)l. the master 
equation for particle (respectively charge) transport, and 
(0 , the equation governing the spin evolution. In the fol- 
lowing, we show that under the approximation of a large 
Rashba length KR t <C 1 , the disorder average of Eq. 
also gives an averaged equation for spin, where the av- 
eraged equations have the same relation to each other, 
as is the case for Eqs. Q and @- In order to elucidate 
the correspondence between the disorder averages of both 
equations, in the following, we first consider a "generic" 
averaging procedure for charge transport, and then ex- 
plore the outcome of the same procedure as applied to 
spin transport. 

The "generic" averaging procedure follows the proce- 
dure used by Klafter and Silbeyi^ which they employed 
to study the connection between the continuous time ran- 
dom walk approach and the exact (generalized) master 
equation. We assume, that the spatially disordered sys- 
tem can be adequately represented using an ordered host 
lattice with very small lattice constant. Then only some 
sites of the host lattice correspond to the original sites 
and are available for transport, whereas all other host 
lattice sites are unavailable in a specific disorder real- 
ization. The transition rates W mim (where m and mi 
correspond to sites of the original [disordered] model) 
can be generalized to transition rates V nin on the host 
lattice, such that V nin vanishes, except when both sites 
are available as defined above. Using this ordered rep- 
resentation of the originally disordered problem, one can 
construct a formal solution and thereafter perform the 
disorder average 

(seeRef.GJ). 

In this way, one obtains a 
generalized master equation (GME) for the (averaged or 
macroscopic) disordered system, which is formally exact. 
Approximations (as, e. g., t he continuous time random 
walk considered in Ref. |l8|) are only needed thereafter, 
in order to obtain explicit expressions for the transition 
rates of the GME and thus be able to calculate solutions 
of the GME. In the following paragraph we will show 
that, provided KR t <C 1, the derivation of the formally 
exact GME for charge transport [starting from Eq. 
can nearly unaltered be applied to spin transport [start- 
ing from Eq. (J2J], too. In this way, any approximations, 
which are applied to the charge transport GME, yield 
immediately a corresponding spin transport GME with- 
out the need to introduce further approximations (except 
KR t <C 1 of course). 

We denote the (disordered) transition rates on the or- 
dered host lattice (with site indices n, etc.) by 



and introduce diagonal elements by 

Kin ^ y l^mmi ^7] 



(Al) 



(A2) 



where S nm is equal to unity or to zero, when the host 
lattice site n coincides or not with an original site m in 
the current disorder realization, respectively. Then, the 
charge transport rate equation can be written as 



dt 



Pn{t) 



(A3) 



These equations can equivalently be expressed as a ma- 
trix equation, such that, if N denotes the number of sites 
of the host lattice, the p n become an N- vector p and cor- 
respondingly V ni n becomes an N x TV-matrix V_. The rate 
equation then takes the simple form dp/dt = V_p, or, in 
Laplace space [s — V_]p — p Q , where the index marks 
the initial conditions. The averaged solution is formally 
obtained as 



(p(s)) = ([ S -Y]- 1 )p () = [s-m^ 



(A4) 



where (...) denotes disorder average. The second equa- 
tion constitutes the definition of the self-energy M, the 
elements of which are the (as yet exact, but only formally 
defined) transition rates of the GME 



Sp n (s) - }^M nin (s)p ni (s) 



PnO- 



(A5) 



Conservation of probability allows to write this equa- 
tion in the usual form 

SPn(s) - PnO = }}M nin (s)p ni (s) - M nm (s)p n (s)] , 

where the sum now excludes the term n\ — n. Note, that 
this equation lives on an ordered (quasi-continuous) lat- 
tice in contrast to Eq. Q). One can see, that the formal 
structure of Eq. Q has survived the averaging, except 
that the transition rates have become frequency depen- 
dent. Thus, even though the microscopic equations 
are Markovian, the macroscopic (averaged) equations can 
be non-Markovian. 

Let us now execute the same procedure for the spin 
density. The quantity p n is a 3-vector, thus p lives in 
the Cartesian product space of 3-vectors and TV-vectors. 



The transition rates are products D nin V ni , 



V„ 



where D nin determines the action in 3-space and V nin 
is identical to the corresponding quantity defined above 
for charge transport. Note, that we do not need a spe- 
cial rule for the "diagonal elements" , because D nn is the 
unity matrix. Thus, the formal solution is 



(p(*)) 



■n 



[s-iwr 1 ^. (a?) 



defining the quantity M . 

We now introduce the assumption, that a product of 
several rotation matrices D obeys the relation 



D„ 



D r , 



(A8) 
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This equation is valid to first order in K. Therefore, the 
condition KR t <C 1, where R t is a relevant length scale 
of the hopping transport, assures the applicability of the 
approximation. The fact, that Eq. I|A8(I is valid only to 
first order in K, gives the reason for the restriction to 
two-site hopping probabilities in this paper. Three-site 
probabilities yield interesting physics 8 ' as, e.g., the spin 
Hall effect, but those effects are of higher order in K. 

Within the approximation l|A8(l . it can be shown, that 
the spin dependence of the transition matrix M_ is given 
by — 



sites M„ 



M(R nin ). Thus, it is convenient to work 



(A9) 



and the rate equation of the averaged system reads 

s Pni s )-PnO = ^[Mmn&Dnm-PmisyMnm&Pnis)}. 

m 

(A10) 

The approximate validity of Eq. I|A9(1 can be seen by 
representing the inverse matrix [s — (VQ] -1 occurring in 
Eq. (|A7|) as a geometric series. In each term of this series, 
the product (|A8() occurs, and is accordingly simplified. 
The disorder average then yields an expression with the 
structure 



(All) 



The series of averaged powers of V_ are equal to the ex- 
pression [s — M] from charge transport, which can sim- 
ilarly be expanded into a geometric series, and the overall 
l)-factor split [as if Eq. (|A8|) is read from right to left] , so 
that each term D nin (M_ )mn of the new geometric series 

becomes (M_ ) nin and we thereby obtain Eq. I|A9(1 . 

Another way to derive the approximation, which uses 
the Zwanzig projection operator method, is the follow- 
ing. Introducing a projection operator P as effecting the 
disorder average 



PA = (A), 



(A12) 



the self-energy can be expressed as (see, e.g., Ref. ITsI) 



n -1 



K = (£) + (sz [s - (i - p)Z\ 6 Z/ ' ( A13 ) 

where 8V_ = V_ — (V_) . Here, only one term, namely [s — 

(1 — P)V}~ 1 ■, has to be expanded in a geometric series. 
The application of approximation (|A8|) again leads to the 
conclusion (|A9() . 

Homogeneity of the averaged system implies, that 
M nin only depends on the difference vector between the 



in wave-vector space. The disorder averaged rate equa- 
tions in Fourier-Laplace space then become the Eqs. © 
and (QJ, where we have replaced the symbol M by W, 
with the understanding that the quantities W mim are 
the transition rates for the disordered system, whereas 
W^(s|q) or W(s|r) are the transition rates for the aver- 
aged — and thus homogeneous — system. 
APPENDIX B: THE RATE EQUATIONS IN THE 
LONG WAVE-LENGTH LIMIT 

The following expressions are written in the right- 
handed orthonormal basis e r = r/r, &Kxr 

= K x r/Kr, 

and e z = K./K. 

The spin rotation matrix for a transition between sites 
m and mi depends only on the distance vector between 
both sites D mim — D(R mim ). The explicit expression 
for this quantity ism 

D(r) = J3 + sin(2Kr)(e r o e z — e z o e r ) 

+ (cos(2A"r) — l){e z o e z + e r o e r ). (Bl) 

First, we show that the Fourier transform l)(q) has 
finite support, specifically, it is restricted to the disc 
|q| < 2K. The (lengthy) calculation is not given here 
in detail, but proceeds in the following way: The (di- 
vergent) Fourier integral of D(r) is made convergent by 
introducing the factor e~ ar with a positive parameter a. 
The corresponding transform _D Q (q) can be calculated 
analytically and consists of rational functions of q and 
a, except for non-analytic factors of [q 2 + a 2 ] -1 / 2 and 
[q 2 + (a±2Ki) 2 ]~ 1 / 2 , which only have cuts and singular- 
ities in the disc q < 2K. The subsequent limit a — > + 
shows that D (q) is a sum of ^-distributions and their 
derivatives with support within this disc. Specifically, 
for q > 2K one obtains lim Q ^ + D a (q) — 0, the zero 
matrix. Thus, 



D(q) = for q > 2K. 



(B2) 



Next, the first and second spatial derivative of D(r) are 
calculated. Recall, that the vector r is two-dimensional. 
Thus, e.g., V o r — e r o e r + &Kxr ^Kxr, where one 
should note, that the z-components are missing. Keeping 



this in mind, we have V o e r 



eKxr/r. Of 



course V o e z = 0, and thus V o e^ xr = V o (e 2 x 
e r ) = — &Kxr e r /r. Replacing the dyadic product by 
a scalar product (corresponding to a tensor contraction) 
immediately gives the additional relations V • e r = 1/r 
and V • exxr = 0. 

One can now easily calculate 
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V o D(r) 

and 

AD(r) = 



2K cos(2Kr)e r o (e r o e z — e z o e r ) H xr o (e^xr ° e z — e 2 o e^ xr ) 

r 

cos(2a>) — 1 

— 2K s'm(2Kr)e r o (e z o e z + e r o e r ) H ejf Xr ° (ejf Xr o e r + e r o ejf Xr ) (B3) 



, , . cos(2A>) sin(2Kr)\ , 

-4K 2 sin(2A> + 2K - - y — — - (e r o e z - e z o e r ) 

r r z J 

, r _ * n „sm.(2Kr)\ . , cos(2A>) - 1 . 
4a cos(2a>) + 2i\ I (e z o e z + e r o e r j — 2 (e r o e r — e/fxr ° ejf> 



(B4) 



r 



The limit r->0 needs a special comment, because one 
formally obtains expressions as, e.g., 

V o D(r)| r= o = 2Ke r o (e r o e z — e z o e r ) 

+ 2Ke Kxr o (e Kxr o e z - e z o ejf X( .), (B5) 

which contain the basis vectors e r and e^- xr , which are 
indeterminate in the limit r = 0. This seeming puzzle is 
solved by the observation, that these basis vectors only 
occur in combinations, which yield the unit matrix I2 in 



two dimensions. 
Thus, one obtains 

a- VoZ)(r)| r=0 = 2K(aoe z e z oa), (B6) 

where a is an arbitrary vector, and 

AD(r)\ r=0 = -4K 2 I 3 -4K 2 e z oe z . (B7) 

Inserting these expressions into Eq. © gives Eq. (JSJ). 
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FIG. 1: Temporal evolution of the components of the total spin polarization for a topologically disordered hopping system from 
numerical simulation: p x (t) (filled squares), p y (t) (filled triangles), p z {t) (filled circles). For comparison the values computed 
from the analytical expressions 11411 and 115L where D and p are replaced by the corresponding time-dependent quantities 
D e g(t) and p e s(t) [obtained numerically from the charge transport in the same ensemble of disordered systems], are also 
displayed: p x (empty squares), p y = 0, p z (empty circles). The statistical error is of the order of the symbol size. The 
parameters are aN~ 1/2 = 3, e = 20, and K = s/N /10. 
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FIG. 2: Temporal evolution of the z-component of the total spin polarization for three values of <x/V _1//2 = 3 (full circles), 4 
(full triangles), 5 (full diamonds). For comparison, the analytic solution for aJV^ 1 ^ 2 — 3 is also given (empty circles). A quite 
high value of the electric field (e = 200) has been chosen in order to numerically reach the oscillatory regime (Note the negative 
value of p z at large times, which is not possible for plain decay.). In contrast to Fig. the time axis is logarithmic. Thus, with 
increasing disorder the spin life-time is strongly enhanced. K = \/~N /10. 
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FIG. 3: Temporal evolution of the z-component of total spin polarization for e = 0. Numerical simulation (full symbols) for 
olM~ x I 2 = 3 (circles), 4 (triangles), 5 (diamonds) and analytical calculations (empty symbols, respectively) are shown. (On 
this scale, the empty diamonds are indistinguishable from the full ones.) K — \/Af /10. 
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FIG. 4: Time dependence of the effective diffusion constant D e g(t) for disorder parameters aH -1 ^ 2 = 3 (circle), 5 (diamond), 
and 7 (square), numerically determined from charge transport in a topologically disordered hopping system. The scale of the 
ordinate is arbitrary. 




FIG. 5: Spin components in the stationary state for a disordered hopping system with spin injection simulated by the boundary 
condition p(x = 0) = e z . The results of numerical simulation p x (square), p y (triangle), p z (circle), compare well with the 
analytical results Eq. 1181 : p x (dotted), p z (dash-dotted). The {/-components of the spin are zero within the statistical error. 
The parameters aM' 1 ^ 2 = 7, e = 0.232, and K — \/Sf /20. The dimensionless length x — |r • e x \/K measures the distance from 
the boundary and the in-plane electrical field is aligned perpendicular to the boundary at x = 0. 
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FIG. 6: The x-component of the spin in the stationary state for three values of ct/V -1 ^ 2 = 3 (circle), 5 (diamond), 7 (square), 
and the analytical result Eq. (1181 (full line). The situation is analogous to Fig. [3] except for a stronger electric field e = 232. 
One can see, that for large e, the magnitude of the spin decreases with increasing disorder, and differs increasingly from the 
analytical results. (The z-component behaves analogously.) A discussion of this behavior is given in the text. 



